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£N| ' We use a lattice Boltzmann method to study pattern formation in chemically reactive binary 

>V 

fluids in the regime where hydrodynamic effects are important. The coupled equations solved by 
the method are a Cahn-Hilliard equation, modified by the inclusion of a reactive source term, and 
the Navier-Stokes equations for conservation of mass and momentum. The coupling is two-fold, 
resulting from the advection of the order-parameter by the velocity field and the effect of fluid 
composition on pressure. We study the the evolution of the system following a critical quench for 
a linear and for a quadratic reaction source term. Comparison is made between the high and low 
viscosity regimes to identify the influence of hydrodynamic flows. In both cases hydrodynamics is 
found to influence the pathways available for domain growth and the eventual steady-states. 

o . 

O , PACS numbers: 47.54.-r, 47.70.Fw, 82.40.Ck 



I. INTRODUCTION 

o 

The process of phase separation in chemically reactive mixtures has been considered by several authors. Glotzer et 
al ~\ and Christensen et al used a modification of the Cahn-Hilliard equation to investigate the effects of a linear 
reaction of the type A <-» B occurring simultaneously with phase separation following an instantaneous quench. In 



o 



o 
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contrast to phase separation alone, domain coarsening was halted at a length-scale dependent on system parameters 



tins 



■ to an effective long-range repulsion 

an effective long-range interaction arises because of an interplay between interactions and steric constraints 



is wn 



resulting in the 'freezing in' of a spatially heterogeneous pattern. It was recognized that the steady-states resulted 
from competition between the demixing_effects of phase separation and the equivalence of the chemical reaction term 

Similar physics is seen in the phase ordering of block copolymers where 

In 

such systems pattern formation is a result of thermodynamic equilibrium. By contrast, in the systems we consider, 
the steady-states are maintained dynamically by the interplay of reaction and diffusion. 

A number of chemically and structurally more complicated systems have been considered, numerically and theoret- 
ically, within the same framework of a modified Cahn-Hilliard equation. These include ternary mixtures 
and systems with orientational order 



Here we investigate the effect of hydrodynamic interactions on phase ordering in a binary fluid mixture with chemical 
reactions using a lattice Boltzmann method. The case of the linear reaction has been considered before by Hou et al 
\l\ by a different numerical method. We duplicate some of their results as a means of testing our approach and then 
consider the quadratic reaction mechanism A + B <-> IB. 



The inclusion of hydrodynamics is known to strongly affect the way in which an unreactive fluid mixture coarsens 
in the aftermath of a quench The growth exponent is found to increase from a = 1/3, for the purely diffusive 

case, to a = 1 or a = 2/3 for the viscous and inertial hydrodynamic regimes respectively. The new pathway for 
growth provided by hydrodynamics is transport of the bulk fluid down a pressure gradient established by variations 



in curvature 
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In two dimensions this minimises curvature by making domains circular, whereupon the effect 



vanishes and further coarsening can only occur by diffusion In addition there is the possibility, investigated by 



Tanaka 
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], that the rapid decrease in interfacial area resulting from the hydrodynamic mechanism may leave the 
bulk phases unequilibrated and subject to a round of secondary phase separations. This suggests that coupling a 
modified Cahn-Hilliard equation to the Navier-Stokes equations for fluid flow may uncover behaviour different to that 
observed for the purely diffusive case. 



Experimental work 



ly Qinus] 
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has shown that a variety of mesoscopic structures can be formed when chemical reactions 



are photo-induced in phase separating polymer mixtures. The effects of two kinds of photo-chemistry have been 

and intramolecular photoisomerisation ^ , ^| . Both give rise 



er mixt 



considered: intermolecular photodimcrisations 
to a long-range inhibition which prevents phase separation proceeding beyond a certain domain size. In the first case 
the inhibition is due to the formation of a network of cross-linked polymer molecules whereas in the second case it 
arises from the differing chemical properties of the two isomers. The similarities in the patterns formed due to phase 
separation arrest in simple fluids and in reactive polymer blends suggest the latter may be approached by considering 
first a small-molecule system. 

The paper is organized as follows. In Section [H] we present a model of a chemically reactive binary fluid which 
couples the processes of reaction and diffusion to flow. We then outline the linear theory of pattern formation in the 
absence of hydrodynamic effects. In Section ITTT1 we construct a lattice Boltzmann scheme which solves the equations 
of motion of Section [H] in the continuum limit. In Sections IIVI and Ivl results are presented for the evolution of both 
high and low viscosity systems after a critical quench for a linear and a quadratic reaction mechanism respectively. 
For the reaction of type A <-> B, comparison is made with the results of 



auadratic read 
, and Q. 
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II. THE MODEL 



A. Equations of motion 



In this section we summarize a model which describes the phase behavior and hydrodynamics of a two-component 
fluid. Labeling the components A and B, we choose a description of the fluid in terms of the following variables: the 
total density, p — pa + Pb', the total momentum, pu, and a compositional order-parameter, <f> = pa — Pb- 

The composition of the fluid evolves according to a modified version of the Cahn-Hilliard equation which includes 
the effects of chemical reaction; advection of the order-parameter by the flow-field, u, and diffusion in response to 
gradients in chemical potential: 



Here Mq is a mobility constant and J, which depends on the reaction rate constants, is the change in 4> per unit time 
due to chemical reactions. The chemical potential of the system, p, is given by the functional derivative of the free 
energy, J 7 , with respect to <fi. 
We choose a free energy 



7 is taken to be greater than zero for stability and the sign of e determines whether the polynomial contribution 
to the free-energy density has one or two minima, and hence whether the fluid is above (e > 0) or below (e < 0) 
its critical temperature. For e < the mixture will separate into two bulk components separated by a narrow, but 
smooth, interface. The gradient-squared term in <f> associates an energy cost with variations in composition and the 
parameter k is related to the surface tension and governs the width of the interface between the two phases. The 
parameter T appears in the isotropic part of the pressure tensor and is related to the degree of incompressibility of 
the fluid ^(J. A suitable choice is T = 1/3. 

We consider two types of reactive source term, J. A linear source 



9*0(x, t) + d a (f>u a = M dlp(x, t) + J(x, i). 



(1) 




(2) 



JM = p(r 2 -r 1 )-0(r 1 + r 2 ), 



(3) 



corresponding to the reversible chemical reaction A «-» B. And a quadratic source 




(4) 



corresponding to the reversible chemical reaction A + B <-> 2B. The constants Fi and I^ are the rates of the forward 
and backward reactions respectively. We note that, for a spatially homogeneous system, the linear mechanism has a 
single stable fixed point whereas the quadratic mechanism has a stable fixed point at <p = (T2 — ^\)p/(J?i + T2) and 
an unstable one at <p = p. Here we consider only cases where ri = Y2 = T. 
The velocity field obeys a Navier-Stokes equation, 

d t pu a + dppu a up = -dpPap + h>dpS a p, (5) 

where P a p is the pressure tensor, v is the viscosity and S a p is the viscous stress tensor. The pressure tensor is derived 
from the free-energy: 

dpP afi = pd a ( S -P\ + <\>d a . (6) 



8p J \ 6<p J 

This provides a further coupling between the evolution of <fi and u in addition to the advection term in ■ From © 
and J2J it follows that 



P, 



a/3 



T P + 0fy/„ - /o - «0V 2 - - (V^) 2 j S al3 + Kd a <j>d <t>, (7) 



where /q = ecf) 2 /2 + 7</» 4 /4 denotes the polynomial contribution to the free-energy density. 
The total mass density of the fluid p is also conserved and obeys 



d t p + dppufi = 0. (8) 



B. Stability and steady-states without hydrodynamics 

Linear stability analysis of the reaction-diffusion equation (JTJ with u = and source term (J3J) shows that only those 
modes, </>(k), with fc; < |k| < k u are unstable, where k u and ki depend on the parameters e, 7 and Mq in equation Q 
\l, ^]. This is in contrast to spinodal decomposition (the case J = 0) where only short- wavelength modes are stable. 
The damping of long-wavelength modes in the reactive case prevents continued growth of domains. Instead, phase 
separation is halted at some length-scale set by the reaction. In addition, there is a threshold value of T, 

T th - M e 2 /Sk, (9) 

above which the reaction is strong enough to completely inhibit phase separation by rendering all linear modes stable. 



We also note, followine 



ing 3 and j], 



that there is an equivalence between this behaviour and phase ordering in a 



system with competing short and long-range interactions. The equivalence can be demonstrated by incorporating the 
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reactive source into an effective free-energy for the system where it appears as a non-local term. Formally, we rewrite 
(1) as 

d t 4>(x, t) + d a (f>u Q = M dl 

where 

and G(x, y) is the Green's function of the Laplace operator. The reaction is then seen to act as an effective long-range, 
intra-species repulsion, with strength governed by T, in contrast to the short-range attraction of like-molecules which 
drives phase separation. 

For the case of equal forward and backward rates, the linearised behaviour of the quadratic source (@J) with reaction 
rate T is the same as for the linear source with reaction rate Tp/2. This can be seen from the linearisation of Q 
with u = and source term Q. Hence, at early times, the linear and quadratic cases segregate in the same way. 
However, after the formation of separate A and B-rich regions, non-linear contributions to the source term become 
important and this is expected to lead to growth of the A-rich phases at the expense of the B-rich ones. This follows 
from the asymmetry of J between the two phases: in the A-rich phase the production of B is limited by the amount 
of B already present. Hence it is limited by the minority component of the fluid in these regions. In the B-rich regions 
the production of A is limited by the amount of the majority phase present. Therefore production of A in the B-rich 
phase is the more rapid process. 



{T-H) 



(10) 



III. A LATTICE BOLTZMANN SCHEME 



The lattice Boltzmann method is a well-established numerical technique for hydrodynamic problems |21| . Initially 
it was a kinetic-theory based method for the simulation of isothermal ideal flows which was introduced to circumvent 
some of the problems which rendered its predecessor, lattice-gas cellular automata, impractical. However, it has since 
been modified and applied to a variety of problems in the simulation of complex fluids. Examples include binary fluids 
[22 . 2^1 24^ liquid-gas systems 2^, 2^, 2^|, liquid-crystals [^J and colloidal suspensions 



A lattice Boltzmann 



scheme for the simulation of two or more species undergoing reaction and diffusion in a moving, viscous solvent was 
formulated by Dawson et al |27| . In comparision, our model incorporates the thermodynamics of the multi-component 
fluid via the Cahn-Hilliard equation. 



FIG. 1: (Color online) The fifteen-velocity face-centered cubic lattice. The dots represent lattice sites and arrows the velocity 
vectors used in the lattice Boltzmann scheme. The lattice spacing, Sx, and time-step, St, are set to unity. 

To simulate the binary fluid model described in Section ^ we utilise the free-energy lattice Boltzmann method of 
Swift et al To this end we define two populations of dynamical variables {/i(x)}" =0 and {<?i(x)}™ =0 on the sites 
of a simple lattice in three-dimensions. On each site the variables /, and gi correspond to a velocity direction e, for 
i = 0, 1, .., n. The dynamical variables are referred to as distribution functions since their moments over the velocity 
set define the macroscopic physical quantities: 

n n n 

P = ^2fi(x,t), pu = ^2eifi(x,t), </> = ^2gi(x,t). (12) 

i=0 i=0 i=0 

The distribution functions on each site are updated in discrete time with a time-step St. The velocities are chosen 
so that eo = and, for all i ^ 0, &i8t lies between two lattice sites. The choice of lattice and velocity set are subject to 
certain restrictions 21] . For this work we used a face-centred cubic lattice in three dimensions with the set of fifteen 
velocities, i = 0, 1, .., 14, illustrated in Figure^ 

The distribution functions fi and gi evolve according to 

/i(x + 5te i; t + St) = ffc t) + A T / [/], (13) 

<?,(x + Ste t ; t + 5t)= . 9l (x; t) + A[ 9 [g] + F u (14) 

where 

AT[/] = -^(/ i (x;i)-/f(x;t)), (15) 

and the relaxation times r/ and r g are free parameters. Equations i|13|) and l|14|) are both lattice equivalents of the 
Bhatnager-Gross-Krook (BGK), or single relaxation time, approximation to the full continuum Boltzmann equation 



We now need to specify the local equilibria functions f^ q and g\ q and the forcing term F\. f^ q and g eq are taken to 
be series expansions in the velocity 



fi 9 = A i + B i e ia u a + Qu 2 + Die ia eipu a up + E" 13 e ia e ll3 , 
g eq = Hi + Jie la u a + KiU 2 + Qie ia eipu a Uf3. 

The coefficients in and (|17|) are chosen so the moments of the equilibrium distributions satisfy 

n n 
n 

^ fi' ! e ta e l = P al 3 + pu a up, 

i=0 

n n 

i=0 i=0 
n 

Yj9T e iaZi(i = DfJ,S af 3 + (j>U a Up. 
i=0 

In addition, the lattice forcing term is chosen to obey 

n n 



i=0 



One possible choice of the coefficients, such that constraints (|18|) - (|22H hold, is given by 



Al ~ 14 = 30 Pac " A °=P- UA i' 



B 7 _ 14 =p/24, fl 1 _ 6 =8S 7 , 
C 7 _i4 = -p/24, C = 16C 7 , Ci_ 6 = 2C 2 , 
£> 7 _ 14 = p/16, D 1= 8D 2 , 



E a(3 - — 

^7-14 - lg 



E 



a/3 
1-6 



8££ 



a/3 



1 



#1-14 = j^Dp, H 



- MFi, 

J 7 _14 = ^/24, J!_ 6 =8J 7 , 

K r -u = -4>/2A, K = 16K r , K X _ 6 = 2K 7 , 

Q 7 _i4 = 0/16, Qi_ 6 = 8Q 7 
F n = 5t l -J, F X - 6 =5t±- A J, F 7 _ 14 =St^J. 



(16) 
(17) 



(18) 
(19) 
(20) 
(21) 



(22) 



(23) 
(24) 
(25) 
(26) 
(27) 

(28) 
(29) 
(30) 

(31) 
(32) 

(33) 
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The constraints (|18fl - <|22[) ensure that, on length and time-scales large compared to the lattice-spacing and time- 
step, the evolution of the moments (|12[) satisfies the partial differential equations set out in Section Til Al To check 
that this is indeed the case the task of reducing the description of the dynamics in terms the distribution functions to 
one in terms of their moments must be addressed. The reduction can be performed by a Chapman-Enzkog expansion 
of equations (|13|l and (|14|l . Since the details of this are essentially no different to those found in we present only 
the result here. 

The zeroth moment of the fi satisfies equation (JSJ) for conservation of mass. The first moment of the fi satisfies 
the Navier-Stokes equation © with v = \ (t/ — 1/2) St and 

S a /3=p(d a up + d/3U a )+£(p,u), (34) 

where the £ denotes unwanted error terms 2^ . 

The first moment of the gi satisfies the reaction-diffusion-advection equation 



d t (j) + d a (j)u a = ujg 



Ddin - dp ( ^d a P af3 



■LJ g d a (u a J) + J[4>], (35) 



where uj g = St(r g — 1/2) and D is defined in (|21l) . Equation l|35|l corresponds to equation Q with mobility Mq = uj g D 
but with two extra terms. The term in gradients of the components of the pressure tensor is present in some other 
free-energy lattice Boltzmann methods and has been shown numerically to be small in comparison to the desired 
terms 23j. The term in gradients of the reactive source can be seen, to first order in <fi, to be a correction to the 
advecting velocity of order u> g T . For our choices of parameters, this contribution is small in the low viscosity regime 
where the advecting flow field is important. 

This completes the specification of our lattice Boltzmann method. Although the model is inherently three- 
dimensional, we consider only its restriction to two dimensions. As an initial condition we choose the total density 
of the fluid p = 1 at each lattice lattice. The near-incompressibility of the fluid ensures that this value remains 
approximately the same at later times. To imitate the conditions following a rapid cooling of a fluid from above to 
below its critical temperature we initialise </>(x; 0) = <5</>(x) where 5(f) is random noise with 50(x) < 0.01 Vx. The 
parameters e and 7 in the free-energy are chosen so that in the unreactive case the fluid phase separates into regions 
where <fi = ±1.0. We choose k — 0.01 to ensure a narrow interface and choose r g = 1.0 and D = 2.0 which fixes the 
diffusion constant Mg = 1.0. We also choose 5t = 5x = 1.0 and the system size L x = L y = 256, throughout. The 
viscosity of the fluid is controlled by varying the relaxation time 77. For Tf — 400, domains in the unreactive fluid 
grow as i 1 / 3 and the fluid can be taken to be in the diffusive regime where hydrodynamic flows are negligible. To 
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FIG. 2: Steady, or near-steady, states for a linear reaction mechanism in the diffusive (ry = 400, left-hand column) and 
hydrodynamic {jf = 5.0, right-hand column) regimes for two different values of T: T = 0.001, top, and T = 0.0001, bottom. 
A-rich regions are shown in white and B-rich ones in black. 

ensure a hydrodynamic growth exponent of 2/3 we choose tj = 5.0. 

IV. LINEAR REACTION 



cnrireij 

n 



We consider first the case of the linear reaction mechanism (jSJ). Figure [3 compares the steady-state structures at 
low and high viscosities for two different values of the rate constant T. 

At high viscosity (Fig. [21 left-hand column) hydrodynamic effects are suppressed and domain evolution is due 
entirely to the interplay between reaction and diffusion. We observe the laminar steady-state morphologies found in 
The average domain size in the steady-state is dependent on the rate of reaction: decreasing T weakens the 
opposition to phase separation posed by the reaction. Hence, for smaller T, the average domain size in the steady-state 
is larger. 

A measure of the characteristic length-scale is given by the inverse interfacial length lj = L x L y /Lj, where Lj is 
defined as the number of lattice sites x such that <f>(x)(f)(x ) < for at least one x a nearest neighbour of x. Figure 
|3fa) shows the time-evolution of li at high viscosity for different values of T. For low rate constants, lj is close to the 
spinodal growth-law, a = 1/3, at early times before asymptoting to a constant value. For a faster reaction, a steady 



length-scale is quickly selected without prior scaling. These results are in qualitative agreement with those of [2| and 
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(b) 



FIG. 3: Log-log plot of the average domain size, measured by the inverse interfacial length, as a function of time at (a) high 
and (b) low viscosity. The reaction rates are F = 0.001 (□); 0.0001 (A) and 10 -5 (o), respectively. The solid straight line 
corresponds to a = 1/3 in (a) and a = 2/3 in (b). 
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At low viscosities hydrodynamics becomes important and pressure driven flows attempt to eliminate local interfacial 
curvature. The effect on phase separation in an unreactive fluid is to increase the rate of domain growth. However, 
the mechanism is such that for a complicated domain structure, conflicting signals will be sent into the bulk phase as 
to the direction in which flow should be established As a result, a complicated flow-field arises. 

For r = 0.001 a labyrinthine pattern of A and B-rich domains, similar to that observed at high viscosity, is formed 
(compare the top-two panels in Fig.[2J). However, the regions of highest interfacial curvature at high viscosity (Fig. [21 
top-left), where two mis-aligned domains meet, can no longer be maintained in the presence of flow. Instead the 



3 



domains curve gently and form concentric arrays of arcs, sometimes enclosing spiral-like structures 

For r = 0.0001, the additional mixing effects of pressure driven flows lead to a greater degree of coarsening before 
the energetic barrier presented by the reaction is reached than in the purely diffusive case (compare the bottom-two 
panels in Fig- • Hence, the domains at low viscosity are larger and more circular compared to the narrow, elongated 
domains found at high viscosity. At later times there is a period of domain creation with smaller domains forming 



3 



inside existing ones. This is analogous to the 'interface-induced secondary quench' described by Tanaka [19j. Rapid 
hydrodynamic growth has resulted in a state where the processes of reaction and diffusion have not yet brought 
about bulk phases with equilibrium values of cf>. Consequently inhomogeneities in the bulk may grow, triggering the 
formation of new domains within an existing one. Figure Elb) shows the evolution of the inverse interfacial length at 
low viscosity for three different reaction rates. For T = 10~ 4 and T = 10~ 5 , the introduction of new interface is clearly 
seen as a peak in lj. At late times a balance is eventually struck between the flow, diffusion and the reactive repulsion. 
The result is the hierarchy of roughly circular interfaces separated by narrow channels shown in Fig. [^bottom- right). 
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FIG. 4: Time evolution of the domain structure for a quadratic reaction mechanism with rate F = 0.0001 for high viscosity 
(r/ = 400, left-hand column) and low viscosity (r/ = 5.0, right-hand column) at times i = 10,000 (top), t = 20,000 (middle) 
and t — 200, 000 (bottom). The ^4-rich regions are shown in white and the B-rich regions in black. 

V. QUADRATIC REACTION 



For the case of the quadratic source term Q , the steady state is an array of £?-rich domains in an A-rich matrix 
where, in both phases, the order-parameter does not take a thermodynamic equilibrium value. The shape of the steady- 
state domains, their density and the evolution to the steady-state depend on the viscosity and the rate of reaction. 
The left-hand column of Figure 01 shows the evolution of the domain structure at high viscosity for T = 0.0001. 
At early times there is a period of domain growth. However, due to the asymmetry of the reaction mechanism, the 
growth of the _B-rich domains is halted at a certain size after which the £?-rich domains are depleted by the reaction 
with smaller domains vanishing completely. As this process occurs some of the larger domains break-up to reduce 
the surface tension. This creates small, circular domains which either vanish or assume a steady size. At some 
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FIG. 5: Time evolution of the domain structure for a quadratic reaction mechanism with rate F = 0.001 for high viscosity 
(t/ = 400, left-hand column) and low viscosity (r/ = 5.0, right-hand column) at times t = 10,000 (top), t — 40,000 (middle) 
and t = 200, 000 (bottom). The A-rich regions are shown in white and the B-rich regions in black. 

point a balance is reached between the depletion of the i?-rich domains and their growth due to B, produced in the 
surrounding A-rich fluid, being supplied by diffusion. 

If the viscosity is decreased the initial rate of growth is higher as shown in the right-hand column of Figure^] As at 
high viscosity the growth of i?-rich domains is eventually reversed by the reaction. However, as the B-rich domains 
are depleted hydrodynamic flows continue to drive coalescence of nearby domains faster than is possible by diffusion 
alone. Consequently, at intermediate times, there are fewer domains and they are less elongated than at high viscosity. 
At later times, gradients of <p within the A-matrix trigger several secondary phase separations, increasing the density 
of domains at the latest times simulated. By t ~ 200, 000 an approximately steady state has been reached. 

For a faster rate of reaction, T — 0.001, 'worm-like' phases are able to persist until the latest time simulated as 
shown in the left-hand column of Figurc[S] This indicates a transition to a reaction-dominated regime where structures 
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10 3 10 4 10 5 10 3 10 4 10 5 



FIG. 6: Log-log plot of the average domain size measured by the inverse interfacial length for (a) V = 0.001 and (b) T = 0.0001 
and values of tj = 400 (high viscosity, A) and 77 = 5.0 (low viscosity, o). The solid straight lines correspond to a = 1/3 and 
a = 2/3. 

with high curvature, which are unfavoured by now and diffusion, can persist for a long time. 

If the viscosity is decreased, as shown on the right of Figure the structures formed are similar to those at high 
viscosity. However, secondary phase separations are triggered within the A-rich phase increasing the number of small 
domains. 

Figure H3 shows the evolution of the inverse interfacial length for different values of T and v. It can be seen 
that when hydrodynamic effects are present the rate of reduction of interfacial length is faster. Moreover, lj is no 
longer monotonically increasing: it attains a maximum due to the creation of more interface during secondary phase 
separation events. 

VI. CONCLUSIONS 

In conclusion, we have extended the lattice Boltzmann method to study the effect of hydrodynamics on structures 
arising in phase-separating reactive mixtures for two simple reaction mechanisms. We have found that hydrodynamic 
flow significantly alters both the way in which the domain structure in these fluids evolves and the eventual steady 
states of the system. The results obtained for the linear reaction are in agreement with previous work on the subject. 
For a quadratic reaction an asymmetric domain structure was obtained, with the inclusion of hydrodynamic effects 
leading to secondary phase separation within majority phase. 

It would be interesting to extend the model to incorporate viscoelastic effects arising from constituent molecules with 
internal microstructure: one can envisage chemical processes which can change the local microstructural elements. 
The technique also makes feasible a study of the interaction of reaction and diffusion with imposed flow- fields. 
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